##### This script runs the linear mixed-effects models reported in the Table OA5.

# The following set of models account for the possible BE and WE of country-round contextual variables
# They also account for the compositional effects of individual-level variables

# Model A: two levels with RE at the COUNTRY-ROUND level
Model_A_WE_BE <- lmer(red_sup ~ employed + unemployed + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + stfeco + lrscale + 
                  mean_LFPR_y0 + dev_mean_LFPR_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 + mean_GDP_y0 + dev_mean_GDP_y0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + liberal + mediter + (1 | cntry_round), 
                  data = ess_4_ld_models, weights = pspwght)
summary(Model_A_WE_BE)
check_collinearity(Model_A_WE_BE)


# Model B: two levels with RE at the COUNTRY level
Model_B_WE_BE <- lmer(red_sup ~ employed + unemployed + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + stfeco + lrscale + 
                  mean_LFPR_y0 + dev_mean_LFPR_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 + mean_GDP_y0 + dev_mean_GDP_y0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + liberal + mediter + (1 | cntry), 
                  data = ess_4_ld_models, weights = pspwght)
summary(Model_B_WE_BE)
check_collinearity(Model_B_WE_BE)


# Model C: three levels with RE at ROUND level (3) and then COUNTRY-ROUND level (2)
Model_C_WE_BE <- lmer(red_sup ~ employed + unemployed + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + stfeco + lrscale + 
                  mean_LFPR_y0 + dev_mean_LFPR_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 + mean_GDP_y0 + dev_mean_GDP_y0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + liberal + mediter + (1 | cntry_round) + (1 | essround), 
                  data = ess_4_ld_models, weights = pspwght)
summary(Model_C_WE_BE)
check_collinearity(Model_C_WE_BE)


# Model D: three levels with RE at COUNTRY level (3) and then COUNTRY-ROUND level (2)
Model_D_WE_BE <- lmer(red_sup ~ employed + unemployed + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + stfeco + lrscale + 
                  mean_LFPR_y0 + dev_mean_LFPR_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 + mean_GDP_y0 + dev_mean_GDP_y0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + liberal + mediter + (1 | cntry_round) + (1 | cntry), 
                  data = ess_4_ld_models, weights = pspwght)
summary(Model_D_WE_BE)
check_collinearity(Model_D_WE_BE)


# Model E: cross-classified model
Model_E_WE_BE <- lmer(red_sup ~ employed + unemployed + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + stfeco + lrscale + 
                  mean_LFPR_y0 + dev_mean_LFPR_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 + mean_GDP_y0 + dev_mean_GDP_y0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + liberal + mediter + (1 | essround) + (1 | cntry), 
                  data = ess_4_ld_models, weights = pspwght)
summary(Model_E_WE_BE)
check_collinearity(Model_E_WE_BE)


# Model F: full model
Model_F_WE_BE <- lmer(red_sup ~ employed + unemployed + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + stfeco + lrscale + 
                  mean_LFPR_y0 + dev_mean_LFPR_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 + mean_GDP_y0 + dev_mean_GDP_y0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + liberal + mediter + (1 | cntry_round) + (1 | essround) + (1 | cntry), 
                  data = ess_4_ld_models, weights = pspwght)
summary(Model_F_WE_BE)
check_collinearity(Model_F_WE_BE)


# Model G: two levels with RE at the YEAR level
Model_G_WE_BE <- lmer(red_sup ~ employed + unemployed + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + stfeco + lrscale + 
                  mean_LFPR_y0 + dev_mean_LFPR_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 + mean_GDP_y0 + dev_mean_GDP_y0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + liberal + mediter + (1 | essround), 
                  data = ess_4_ld_models, weights = pspwght)
summary(Model_G_WE_BE)
check_collinearity(Model_G_WE_BE)


# Plotting the model results
tab_model(Model_A_WE_BE, Model_B_WE_BE, Model_C_WE_BE, Model_D_WE_BE, Model_E_WE_BE, Model_F_WE_BE, Model_G_WE_BE, digits = 3, digits.p = 4, digits.re = 4)



#################################################################################
# Storing and exporting results of the seven models into a Table OA5 of the paper #
#################################################################################

# Creating an empty data.frame into which the results will be stored
Models_WE_BE <- data.frame(c("Intercept", "", "employed", "", "unemployed", "", 
                               "hinc_2_cop", "", "hinc_3_dif", "", "hinc_4_vdif", "", 
                               "female", "", "age", "", "isced2", "", "isced3", "", 
                               "isced4", "", "isced56", "", "stfeco", "", "lrscale", "", 
                               "mean_LFPR_y0", "", "dev_mean_LFPR_y0", "", "mean_SOCEXP_y0", "", "dev_mean_SOCEXP_y0", "", 
                                "mean_GDP_y0", "", "dev_mean_GDP_y0", "", "mean_gini_disp_y0", "", "dev_mean_gini_disp_y0", "", 
                               "liberal", "", "mediter", "", "Var_cntry", "Var_essround", 
                               "Var_cntry_round", "Var_individual", "ICC_cntry", "ICC_essround", "ICC_cntry_round"), matrix(data = rep(NA, 55*7), nrow = 55))
colnames(Models_WE_BE) <- c("Parameter", "Model A", "Model B", "Model C", "Model D", "Model E", "Model F", "Model G")

# Using the programmed function to create Table OA5 of the article
Models_WE_BE <- assign_fe_results(model_name = Model_A_WE_BE, assignment_object = Models_WE_BE, column_number = 2)
Models_WE_BE <- assign_fe_results(model_name = Model_B_WE_BE, assignment_object = Models_WE_BE, column_number = 3)
Models_WE_BE <- assign_fe_results(model_name = Model_C_WE_BE, assignment_object = Models_WE_BE, column_number = 4)
Models_WE_BE <- assign_fe_results(model_name = Model_D_WE_BE, assignment_object = Models_WE_BE, column_number = 5)
Models_WE_BE <- assign_fe_results(model_name = Model_E_WE_BE, assignment_object = Models_WE_BE, column_number = 6)
Models_WE_BE <- assign_fe_results(model_name = Model_F_WE_BE, assignment_object = Models_WE_BE, column_number = 7)
Models_WE_BE <- assign_fe_results(model_name = Model_G_WE_BE, assignment_object = Models_WE_BE, column_number = 8)
View(Models_WE_BE)

### Manual assignment of variance components and ICCs
# Model A: RE at cntry_round
Models_WE_BE[Models_WE_BE$Parameter == "Var_cntry_round", "Model A"] <- as.data.frame(VarCorr(Model_A_WE_BE))[, "vcov"][1]
Models_WE_BE[Models_WE_BE$Parameter == "Var_individual", "Model A"] <- as.data.frame(VarCorr(Model_A_WE_BE))[, "vcov"][2]
icc(Model_A_WE_BE, by_group = TRUE)
Models_WE_BE[Models_WE_BE$Parameter == "ICC_cntry_round", "Model A"] <- (as.data.frame(VarCorr(Model_A_WE_BE))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_A_WE_BE))[, "vcov"]))*100

# Model B: RE at cntry
Models_WE_BE[Models_WE_BE$Parameter == "Var_cntry", "Model B"] <- as.data.frame(VarCorr(Model_B_WE_BE))[, "vcov"][1]
Models_WE_BE[Models_WE_BE$Parameter == "Var_individual", "Model B"] <- as.data.frame(VarCorr(Model_B_WE_BE))[, "vcov"][2]
icc(Model_B_WE_BE, by_group = TRUE)
Models_WE_BE[Models_WE_BE$Parameter == "ICC_cntry", "Model B"] <- (as.data.frame(VarCorr(Model_B_WE_BE))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_B_WE_BE))[, "vcov"]))*100

# Model C: RE at essround and cntry_round
Models_WE_BE[Models_WE_BE$Parameter == "Var_essround", "Model C"] <- as.data.frame(VarCorr(Model_C_WE_BE))[, "vcov"][2]
Models_WE_BE[Models_WE_BE$Parameter == "Var_cntry_round", "Model C"] <- as.data.frame(VarCorr(Model_C_WE_BE))[, "vcov"][1]
Models_WE_BE[Models_WE_BE$Parameter == "Var_individual", "Model C"] <- as.data.frame(VarCorr(Model_C_WE_BE))[, "vcov"][3]
icc(Model_C_WE_BE, by_group = TRUE)
Models_WE_BE[Models_WE_BE$Parameter == "ICC_essround", "Model C"] <- (as.data.frame(VarCorr(Model_C_WE_BE))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_C_WE_BE))[, "vcov"]))*100
Models_WE_BE[Models_WE_BE$Parameter == "ICC_cntry_round", "Model C"] <- (as.data.frame(VarCorr(Model_C_WE_BE))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_C_WE_BE))[, "vcov"]))*100

# Model D: three levels with RE at COUNTRY level (3) and then COUNTRY-ROUND level (2)
Models_WE_BE[Models_WE_BE$Parameter == "Var_cntry", "Model D"] <- as.data.frame(VarCorr(Model_D_WE_BE))[, "vcov"][2]
Models_WE_BE[Models_WE_BE$Parameter == "Var_cntry_round", "Model D"] <- as.data.frame(VarCorr(Model_D_WE_BE))[, "vcov"][1]
Models_WE_BE[Models_WE_BE$Parameter == "Var_individual", "Model D"] <- as.data.frame(VarCorr(Model_D_WE_BE))[, "vcov"][3]
icc(Model_D_WE_BE, by_group = TRUE)
Models_WE_BE[Models_WE_BE$Parameter == "ICC_cntry", "Model D"] <- (as.data.frame(VarCorr(Model_D_WE_BE))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_D_WE_BE))[, "vcov"]))*100
Models_WE_BE[Models_WE_BE$Parameter == "ICC_cntry_round", "Model D"] <- (as.data.frame(VarCorr(Model_D_WE_BE))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_D_WE_BE))[, "vcov"]))*100

# Model E: RE at country and ess_round
Models_WE_BE[Models_WE_BE$Parameter == "Var_cntry", "Model E"] <- as.data.frame(VarCorr(Model_E_WE_BE))[, "vcov"][1]
Models_WE_BE[Models_WE_BE$Parameter == "Var_essround", "Model E"] <- as.data.frame(VarCorr(Model_E_WE_BE))[, "vcov"][2]
Models_WE_BE[Models_WE_BE$Parameter == "Var_individual", "Model E"] <- as.data.frame(VarCorr(Model_E_WE_BE))[, "vcov"][3]
icc(Model_E_WE_BE, by_group = TRUE)
Models_WE_BE[Models_WE_BE$Parameter == "ICC_cntry", "Model E"] <- (as.data.frame(VarCorr(Model_E_WE_BE))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_E_WE_BE))[, "vcov"]))*100
Models_WE_BE[Models_WE_BE$Parameter == "ICC_essround", "Model E"] <- (as.data.frame(VarCorr(Model_E_WE_BE))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_E_WE_BE))[, "vcov"]))*100

# Model F: RE at all three effects
Models_WE_BE[Models_WE_BE$Parameter == "Var_cntry", "Model F"] <- as.data.frame(VarCorr(Model_F_WE_BE))[, "vcov"][2]
Models_WE_BE[Models_WE_BE$Parameter == "Var_essround", "Model F"] <- as.data.frame(VarCorr(Model_F_WE_BE))[, "vcov"][3]
Models_WE_BE[Models_WE_BE$Parameter == "Var_cntry_round", "Model F"] <- as.data.frame(VarCorr(Model_F_WE_BE))[, "vcov"][1]
Models_WE_BE[Models_WE_BE$Parameter == "Var_individual", "Model F"] <- as.data.frame(VarCorr(Model_F_WE_BE))[, "vcov"][4]
icc(Model_F_WE_BE, by_group = TRUE)
Models_WE_BE[Models_WE_BE$Parameter == "ICC_cntry", "Model F"] <- (as.data.frame(VarCorr(Model_F_WE_BE))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_F_WE_BE))[, "vcov"]))*100
Models_WE_BE[Models_WE_BE$Parameter == "ICC_essround", "Model F"] <- (as.data.frame(VarCorr(Model_F_WE_BE))[, "vcov"][3]/sum(as.data.frame(VarCorr(Model_F_WE_BE))[, "vcov"]))*100
Models_WE_BE[Models_WE_BE$Parameter == "ICC_cntry_round", "Model F"] <- (as.data.frame(VarCorr(Model_F_WE_BE))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_F_WE_BE))[, "vcov"]))*100

# Model G: RE at essround
Models_WE_BE[Models_WE_BE$Parameter == "Var_essround", "Model G"] <- as.data.frame(VarCorr(Model_G_WE_BE))[, "vcov"][1]
Models_WE_BE[Models_WE_BE$Parameter == "Var_individual", "Model G"] <- as.data.frame(VarCorr(Model_G_WE_BE))[, "vcov"][2]
icc(Model_G_WE_BE, by_group = TRUE)
Models_WE_BE[Models_WE_BE$Parameter == "ICC_essround", "Model G"] <- (as.data.frame(VarCorr(Model_G_WE_BE))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_G_WE_BE))[, "vcov"]))*100

# Exporting Table OA5 to a csv file
write.csv(x = Models_WE_BE, file = "table_OA5_Models_BEs_WEs.csv")



### Collecting the coefficient estimates for the LFPR BE contextual variable and storing them to a data.frame
LFPR_BE <- data.frame(c("Model A", "Model B", "Model C", "Model D", "Model E", "Model F", "Model G"),
                      matrix(data = rep(NA, 7*3), nrow = 7))
colnames(LFPR_BE) <- c("Model type", "Estimate", "Std. Error", "df")

LFPR_BE[LFPR_BE$`Model type` == "Model A", 2:4] <- coef(summary(Model_A_WE_BE))[15, c("Estimate", "Std. Error", "df")]
LFPR_BE[LFPR_BE$`Model type` == "Model B", 2:4] <- coef(summary(Model_B_WE_BE))[15, c("Estimate", "Std. Error", "df")]
LFPR_BE[LFPR_BE$`Model type` == "Model C", 2:4] <- coef(summary(Model_C_WE_BE))[15, c("Estimate", "Std. Error", "df")]
LFPR_BE[LFPR_BE$`Model type` == "Model D", 2:4] <- coef(summary(Model_D_WE_BE))[15, c("Estimate", "Std. Error", "df")]
LFPR_BE[LFPR_BE$`Model type` == "Model E", 2:4] <- coef(summary(Model_E_WE_BE))[15, c("Estimate", "Std. Error", "df")]
LFPR_BE[LFPR_BE$`Model type` == "Model F", 2:4] <- coef(summary(Model_F_WE_BE))[15, c("Estimate", "Std. Error", "df")]
LFPR_BE[LFPR_BE$`Model type` == "Model G", 2:4] <- coef(summary(Model_G_WE_BE))[15, c("Estimate", "Std. Error", "df")]

# Computing the 95% confidence interval (based on student?s t-distribution)
LFPR_BE$lower_bound <- LFPR_BE$Estimate - LFPR_BE$`Std. Error`*qt(p = 0.975, df = LFPR_BE$df)
LFPR_BE$upper_bound <- LFPR_BE$Estimate + LFPR_BE$`Std. Error`*qt(p = 0.975, df = LFPR_BE$df)


### Collecting the coefficient estimates for the LFPR WE contextual variable and storing them to a data.frame
LFPR_WE <- data.frame(c("Model A", "Model B", "Model C", "Model D", "Model E", "Model F", "Model G"),
                      matrix(data = rep(NA, 7*3), nrow = 7))
colnames(LFPR_WE) <- c("Model type", "Estimate", "Std. Error", "df")

LFPR_WE[LFPR_WE$`Model type` == "Model A", 2:4] <- coef(summary(Model_A_WE_BE))[16, c("Estimate", "Std. Error", "df")]
LFPR_WE[LFPR_WE$`Model type` == "Model B", 2:4] <- coef(summary(Model_B_WE_BE))[16, c("Estimate", "Std. Error", "df")]
LFPR_WE[LFPR_WE$`Model type` == "Model C", 2:4] <- coef(summary(Model_C_WE_BE))[16, c("Estimate", "Std. Error", "df")]
LFPR_WE[LFPR_WE$`Model type` == "Model D", 2:4] <- coef(summary(Model_D_WE_BE))[16, c("Estimate", "Std. Error", "df")]
LFPR_WE[LFPR_WE$`Model type` == "Model E", 2:4] <- coef(summary(Model_E_WE_BE))[16, c("Estimate", "Std. Error", "df")]
LFPR_WE[LFPR_WE$`Model type` == "Model F", 2:4] <- coef(summary(Model_F_WE_BE))[16, c("Estimate", "Std. Error", "df")]
LFPR_WE[LFPR_WE$`Model type` == "Model G", 2:4] <- coef(summary(Model_G_WE_BE))[16, c("Estimate", "Std. Error", "df")]

# Computing the 95% confidence interval (based on student?s t-distribution)
LFPR_WE$lower_bound <- LFPR_WE$Estimate - LFPR_WE$`Std. Error`*qt(p = 0.975, df = LFPR_WE$df)
LFPR_WE$upper_bound <- LFPR_WE$Estimate + LFPR_WE$`Std. Error`*qt(p = 0.975, df = LFPR_WE$df)



# Deleting model results from the workspace (as they are no longer needed)
rm(Model_A_WE_BE, Model_B_WE_BE, Model_C_WE_BE, Model_D_WE_BE, Model_E_WE_BE, Model_F_WE_BE, Model_G_WE_BE)

